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This paper describes a method to do ab initio molecular dynamics in electronically excited 
systems within the random phase approximation (RPA). Using a dynamical variational treatment 
of the RPA frequency, which corresponds to the electronic excitation energy of the system, we derive 
coupled equations of motion for the RPA amplitudes, the single particle orbitals, and the nuclear 
coordinates. These equations scale linearly with basis size and can be implemented with only a 
' single holonomic constraint. Test calculations on a model two level system give exact agreement 

I with analytical results. Furthermore, we examined the computational efficiency of the method by 

^\ . modeling the excited state dynamics of a one-dimensional polyene lattice. Our results indicate that 

T-H ' the present method offers a considerable decrease in computational effort over a straight-forward 

f— I , configuration interaction (singles) plus gradient calculation performed at each nuclear configuration. 

^ ■ Pacs: 71.15-m, 71.15.Pd. 31.50.-l-w 
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• I. INTRODUCTION 

I In recent years, there has been considerable advances towards fully ab initio molecular dynamics simulations using 

I : inter- and intramolecular potentials derived from quantum many-body interactions. This activity has been largely 

J — spurred by the introduction of hybrid quantum/classical dynamics, treatments, such as the density functional the- 

PsJ , ory/molecular dynamics methods introduced by Car and ParrinellaJp (CPMD). In this approach, the quantum elec- 

' tronic orbitals, {0^(1}}, are coupled to the classical nuclear degrees of freedom, R„, through a fictitious Lagrangian, 

, which along with an appropriate set of Lagrange multipliers, (Ay ), leads to a set of equations of motion 

6; m„r„ = -Vr£;[{<^j,r„], (i) 

^ where E[{4>i}, R„] is the Hohenberg and Kohn energy functional!. This approach has proven to be highly successful 
Q in simulating chemical dynamics in the condensed phase. 

^ However, there are a number of distinct disadvantages to the method. First, the equations of motion involve both 

ILJ fast (the orbitals) and slow (the nuclei) degrees of freedom. Thus, the computational time-step must be rather small. 
. !^ \ Secondly, imposing the holonomic constraint of orbital orthogonality forces a single iteration of the approach to scale 
' as ~ M X where TV is the number of particles and M the size of the electronic basis. Finally, the method can only 
5h I treat systems in the electronic ground state. 
. . . 1 Several groups have attemptedJ:o overcome each of these difficulties. The time-step problem has been attacked 
via multiple time-scale methodsErQ and by using the electronic density instead of the Koha-Shaim orbitals as the 
dynamical variableH. Several algorithms have been introduced to tackle the scaling problemS. In these cases, scaling 
more favorable than is achieved only at the cost of several assumptions and approximations. 

Recently, Alavi, Kohanoff, Parrinello, aad-Jrenkel have developed a density functional based method for treating 
systems with finite temperature electronajM. This approach has a number of attractive features, foremost being 
the inclusion of fractional occupation of the electronic orbitals. Dynamics in this scheme are isothermal rather than 
adiabatic, allowing for incoherent transitions between electronically excited states. While the approach works well for 
small bandgap systems (such as metals or dense hydrogen), this approach is not useful for systems with larger band 
gaps {AE ^ kT) or photo-excited systems. 

In this paper, we present an alternative formulation of the Car-Parrinello method for treating electronically excited 
systems. We begin by introducing the general computational algorithm for the dynamical optimization of excited 
electronic states in the presence of "classically propagating" ions. This general theoretical scheme is developed without 
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specific reference to the description of the electronic excited state. The principal results of this section are a set of 
dynamical equations of motion for the electronic amplitudes and nuclear coordinates in a generalized phase space. 
We then derive specific equations of motion for the electronic amplitudes and nuclear coordinates under the Random 
Phase Approximation (RPA) for the electrons. Finally, we demonstrate the salient features of the approach through 
a simplified SU(2) model of interacting electrons and via a realistic model of conjugated one-dimensional polymer 
lattices. 

II. CLASSICAL EQUATIONS OF MOTION FOR THE DYNAMICAL OPTIMIZATION OF EXCITED 

QUANTUM STATES 

We begin by developing a general dynamical theory for molecular dynamics in excited electronic states. We first 
define the iV-electron Fock space as being parameterized by a set of dynamical variables A(r) which include the 
nuclear positions and any other dynamical variables and ensemble constraints placed upon the system (e.g. pressure, 
temperature, volume etc.): 



^= ||ni,n2,...(A)),^n, -7v| (2) 

We describe the ground state wave function as an expansion on the basis in T: 

|*„(A)) = a(A)|ni...n,...(A)) , (vI/„(A)|vI/„(A)) = 1 (3) 

|5'o(A)) is not a Slater determinant in general case. 

We next define the mapping operators, Pj^, which act on the ground state to produce a new state, |i^(A)), which is 
orthogonal to |vI'o(A)) in the subspace J-+ of the Fock space: 

|KA))=pt|^^(A)), (4) 
P.|«'o(A)) =0, (5) 

(KA)|*o(A)) = (vI/„(A)|P,|vI/„(A)) = , (KA)IKA)) = 1 (6) 
If t|he 1^*0 (A)) is a true ground state of the iV-particle system, the excited states should lie in the orthogonal subspace 

Variational determinations of the excited states are thus searches in the orthogonal subspace for the optimal mapping 
operator and for the optimal basis of the Fock space \n\n2...ni.}i such that the energy functional 

- (KA)|i/|KA)), 
= (*.(A)|i/|*.(A)) 

+ (MA)|if|KA))-(*o(A)|7?|vI/„(A))), 

= (*„(A)|F|*„(A)) + (vI/„(A)| [P., [H,pt]] |vI/„(A)), 

= -Bo+w^, (7) 

is at a variational minimum. We write the excitation energy, Ej^, as a functional of both the ground and excited state 
and introduce explicitly via the double commutator term the energy gap, Wi/, between the ground and the excited 
state. 

According to the calculus of variations, we can derive equations of motion which dynamically optimize the excited 
state system. We first define formally a set of "velocities" of the states in the Fock space 

\v)^Uy). \^o)^^\^o). (8) 

where r is not the true physical time but rather a parameter which characterizes the evolution of the amplitudes 
through the Hilbert space subject to the orthogonalization constraints given above. The "classical" Lagrangian for 
our system is given by 
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-£.[|0),k),A] 

-A((vi/„|vi/,)-i)-r((H^)-i), (9) 

where A and F are Lagrange multipliers arising from the holonomic ortho-normalization constraint for the ground and 
excited states and which allow and |^'o) to be treated as independent variables. The masses fi and /i' arCpEctitious 
masses and have no true physical meaning. This Lagrangian leads to the classical Euler-Lagrange equationscj for the 
wave functions: 

' 0, (10) 



A^_^ = 0. (11) 

The Euler-Lagrange equation dictates the following equations of motion for quantum dynamical optimization in the 
Hilbert space: 

MV) = -^ + r|.), 

mX^J-^. (12) 
This dynamical view is the starting point for the remainder of our calculations. 

III. RANDOM PHASE APPROXIMATION FOR ELECTRONIC MOTIONS 

For a molecular system with N electrons and n nuclei, the full A^-body Hamiltonian for a fixed set of nuclear 
positions {Ri/}, is given as (atomic units are used throughout) 

= E^^'^ + Eur^' (14) 



ri - r^i 



(i) . 

where ho is the one-body electron operator which includes the electron kinetic energy and the interaction between 
the electron and thej-iiuclear cores. We define the single-particle Hamiltonian„ i7i , through either the Hartree-Fock 
(HF) approximationtj or via the density functional theory of Kohn and ShamQ. In cither case, Hi is a functional of 
the single particle density 

occ 

p(12) = ^0*(lM.(2). (15) 

i 

The single particle states, are solutions of a nonlinear Schrodinger equation, 

77i[p,R]|0,) =i?,[R]|0,), (16) 

which must be determined at the start of the calculation. Here Ei [R] is the orbital energy. The HF spin-orbital wave 
functions and energies differ from the KS ones in the way the terms in the single-particle Hamiltonian are treated. 
We will denote the Slater determinant wave-function constructed from occupied single-particle HF (of KS) orbitals 
as \HF) and note that this state represents either the Hartree-Fock ground state or the Kohn-Sham ground state. 
Having obtained a basis of single-particle states, we can write the two-body Hamiltonian in second quantized form: 
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H = Ehf "-iaO-i-^ ' 

ijkl era' 

where Ei is the HF single-particle eigciicncrgy, 

{ij\v\kl) = j dld2^*{l)<l>*{2)v{12)M^)M^), (18) 

is the matrix element of Coulomb electron-electron interactions, : • • • : denotes a normal ordered product and Ehf is 
the HF ground state energy. The operators al^ {aia) create (annihilate) single electrons in the {icr} HF spin-orbital 
state, where a = ±\ is the spin index of the electrons. 

It is convenient at this point to move to a particle/hole representation and introduce the particle/hole quasi-particle 
operator, 

ai^\HF)=Q, (19) 

through its action on the HF vacuum, where i = p,h represents particles (occupied states above the Fermi energy) 
and holes (vacancies below the Fermi energy). These quasi-particles operators are related to the fermion operators 
written above via 

al = 9{Ei - Ef)al + e{Ef - Ei)ai^, (20) 

a,, = e{E, - Ef )a,„ + e{Ef - E,)al, (21) 

where 6{x) is the Heaviside step function and Ef is the Fermi energy. 

We can expand H in terms of the particle/hole operators and neglect terms with an odd number of particle or hole 
terms (i.e. anharmonic interaction terms of the form a^aaa and a^a^a^a), 

H = Ehf + ^ Epal^apa + ^ Ehahaal^ 

prr ha 



Plh2h3P4 <t<t' 
pih^Pihi era' 

-\- anharmonic terms, (22) 

where p and h refer to particle and hole single particle states. 

Note that each "harmonic" term in our expansion of H involves an even number of particle and hole operators. 
We can form a new set of exciton operators by taking bilinear combinations of the particle/hole operators with total 
angular momentum S and projection mg, 

cri(T2 

Aph{Sms) = y^^{\(yi^a2\Sms)aha2<^hai, (23) 

(Tl(T2 

where {siaiS2(J2\Sms) is a Clebsch-Gordon coefficient. Our notation is such that a denotes the action of time reversal 
operator on the hole state, 
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«L = (-)^-''^«U- (24) 

The {Ap^ (57715), ylp/i(S'ms)} serve to create (or destroy) particle/hole exciton with spin S and projection m^. 
Writing just the "harmonic" part of H in terms of exciton operators, 



per her 
PlP^h^hi Srris 



2 

PlP^hzhi Snis 

2 

pih2h3P4 Srris S' m'^ 

- iPih2\v\p3hi)YKihA^^'^'^^P3h2{Sm,), (25) 

Plh2P3hi Snis 

where Apfi{Sms) = {~)^'^"^'' Aph{S — m) is the time reversal exciton annihilation operator. This defines our "quasi- 
harmonic" exciton Haniiltonian. 

The exact commutation relation between the exciton operators is 

[ Ap^,iSt,),Al,^,iS'fi')] 



(Ti(T2 <y\'(y2' 



We apply the Wick theorem tp-|the right hand side of Eq. and then neglect the normal ordering contribution to 
this relation. This is the RPAllJ and in this approximation the commutation relation between the exciton operators 
is bosonic 

[V(5/i),4,^,(5V)]«<5ss"^ (27) 
We then perform a canonical Bogolubov transformationO for both singlet (5 — 0) and triplet (5 = 1) exciton bosons 

gl (00) = E KhKhioo) - V(oo), 

ph 

Ql(lm,) = E KhKhi^^^) - y;M^^)^ (28) 

ph 

where X^/^ and are the RPA amplitudes which are yet to be determined. These operators create and destroy 
harmonic excitations in the electronic degrees of freedom. Physically, the singlet RPA excitations correspond to 
"breathing" modes of the electron cloud and triple RPA excitations correspond to dipolar oscillations of the electronic 
system about the nuclear frame. In what follows, we will focus our attention upon the triplet excitations, as these are 
the states created upon photo-excitation of a singlet ground state. Equations of motion for singlet RPA excitations 
can be obtained as well. 

The commutation relations for the RPA excitation creation operators goes as follows: 

[Qiii, gjx'i'] = E '^pf^ '^p'^ ~ ^p'^ ^p'^ 

ph 

= Su'S^^>. (29) 



By defining the RPA vacuum as 



and the RPA excited state as 



Q;.,|RPA)=0, (30) 
qI,, |RPA) = (31) 
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the orthogonality condition of Eq. ^ is fulfilled automatically for single RPA excitations 

(RPA|Q^,Qt,^,|RPA) = (RPA|[Q^„Qt,,,]|RPA) 

= 5^ti'5ii'. (32) 

The RPA excitations are true oscillations of electronic systems in the sense that the Heisenberg equations of motion 
for the excitation creation operators is precisely as we expect for a harmonic system, 

[H2,QU^QI,u:,. (33) 

where oji is the RPA frequency. From this wc can derive the RPA amplitudes, X*^ and 1^^, and RPA frequencies oji. 
After a bit of algebra we arrive at the equations: 

Epux;^ + Y,(pp'\v\h'h)Y;,^, - Y,{ph'Wh)xi,^, 

p'h' h'p' 

= ^rX;H, (34) 

E,^,YKph) + Y,{pp'\v\h'h)Xl,^, - Y,{ph'Wh)Y;,^, 

p'h' h'p' 

= -c^.r;,. (35) 



In matrix form: 



A Q\(\X)\_ ( \X) 



(36) 



The diagonal blocks, A, are determined by Coulombic interactions between single particle/hole excitations. 

Aphp'h' = Eph6pp>dhh' ~ {ph'\v\p'h). (37) 

However, the off diagonal terms, 

Bphp'h' = {pp'\v\h'h), (38) 

are related fluctuations about the HF vacuum and thus give rise to a net polarization of the HF vacuum. Thus, the 
RPA vacuum differs from the HF vacuum by the addition of these polarization fluctuations. Lastly, we note, that if we 
take F ^ or assume that B « and neglect the fluctuation terms, the RPA equations reduce to the Tamm-Dankoff 
or configuration interaction equations with single excitations. . .. . 

The matrix RPA method can be also derived as the small ampl itud e limit of the time dependent HF equationESo. 
We emphasize that the time-dependent density functional theoryllU'ta is equivalent to matrix RPA if we also assume 
harmonic oscillation of the time-dependent density p{t) around equilibrium. The standard practical realization of 
the time-dependent density functional theory is the time-dependent local density approximationta where nonlocal 
exchange interaction is approximated by local density dependent interaction in adiabatic limit. In this respect the 
matrix RPA method we use is advantageous because this method can treat nonlocal Fock interactions in natural way. 

We can transform this into an eigenvalue problem by multiplying the right hand side of this equation by the Pauli 
spin matrix, ctz, 

Ml^-) = WCT^I*), (39) 

where 

l*> = ( \y] ) > (40) 

is a two component spinor state. The eigenvector solution to this equation can be obtained by transforming the RPA 
matrix equation 

crjM|*) = tj|*). (41) 
Thus, 1^) is the right-handed eigenvector of the antisymmetric matrix 
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-B -A 

The RPA eigenstates are subject to the orthogonahty condition 



^"M = ( ] . (42) 



'^^ph ^ph - Yph Y^h - ^ij- (43) 

ph 



Thus, we define ^ as the right-handed eigenvector 



and as the left-handed eigenvector 



l*f) = ( ) , (44) 



l*f>-( l^^l )=<^.|*f>- (45) 



The inner product relation between the left and right-handed eigenvectors produces the desired orthogonality relation 
for the RPA states. 

(*.fl*f> -^f = (46) 

These states are determined at the start of the calculation in a basis of single particle/hole states which spans the 
Fock space. 

The RPA frequency, uji, is determined by taking the expectation value of M between the left and right-handed 
eigenvectors. 

c.,;[X,r,(/.,R] = (1'f|M|vI/f) 

A ^ \( Xk 



A. Equations of motion for the excited states molecular dynamics in the RPA 

The energy functional for the RPA equations is 

E,S) = (RPA|Q^,i7QtjRPA), 
= (RPA|iJ|RPA) 



(RPA| 



ff,g^ll I RPA), 



where Eo is the RPA ground state energy and 

uj, = (RPA| [g^„ \ii,Qi^ II^PA), 



(48) 



(49) 



is the RPA frequency. In the previous section the RPA equation is derived via the Heisenberg equation (Eqjs^) of 
motion for the excitation creation operators. The Heisenberg equation is physically the same as the requirement of 
a variational minimum for the RPA excited state. To demonstrate this, we multiply from the left with an arbitrary 
variation of the RPA excited state: 



(RPA|5g^, \h,qI^ |RPA) =cc;,(RPA|,5g^,g^|RPA) 



(50) 



We can rewrite this equation as a variational of the double commutator in Eq. E^, because (RPA|gj^j = 
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^ {(RPA| [q^., [h,qI]] |RPA) -^,((RPA|Q^,Qj;jRPA) - 1)} = (51) 

where coi plays a role of the Lagrange multiplier and insures the normalization of the RPA excited states. 

Neglecting vacuum polarization effects, i.e. approximating |RPA) by the HF Slater determiuant \Qhf) in Eq. |5l| 
we obtain the so-called equation of motion method which was-jOpginally proposed by D.J. RowetJ and widely used in 
the studies of molecular excited states in quantumjfihcmistryEHrEd. The same matrix RPiL-pouation can be derived by 
several methods, e.g. by Green function methodsEJa, linear response function mcthodLJ-L3 Each of these methods 
has the strong points in specific aspects but is deficient in other respects. As we demonstrated by Eq.^ our method is 
equivalent to the Rayleigh - Ritz variational principle for the RPA excitation amplitudes and the use of the variational 
approach is pivotal to formulation classical equations of motion for the dynamical optimization of the excited states. 

Finally, we assume that polarization of the ground state vacuum due to particle-hole interaction is weak, then the 
RPA vacuum is approximately the same as the HF or KS vacuum and we can write the RPA ground state energy as 
Eo » Ehf- 

Following these considerations, we can write the Lagrangian 



i k ph 

k ph 



- EHF[^,Ti]-LJ,[cj),X,Y,Ii] 

+ ^Tki{{^k\^i)~S,,) (52) 



kl 

where the constraint 



5ZrM((*fc I *i 



kl 



— E I E -^ph-^ph ~ '^pvXph 
kl \ ph 

(53) 

insures the orthogonality of the RPA excited states. Solving the Euler-Lagrange equations for the single particle 
amplitudes, ^fc(l), 

J/0i(l) = 



(50* (1) <50*(1) 
+ 5^A,,0,(1) (54) 
j 

The equations for the single particle orbitals are similar to the CP equations written above (Eq. ^) and include a term 
due to the electronic excitations. In the case of a collective excitation, in which the RPA excitation is delocalized over 
a number of single particle states (or in other worlds, the oscillator strength of the RPA excitations is distributed 
more or less uniformly over a number of particle-hole excitations) the variation of the energy of collective excitation 
Lu with an infinitesimal change in one of the single particle orbitals will be very small. Because of this, the RPA 
excitation variables will evolve on a slower time-scale than the single particle variables and we can invoke an adiabatic 
separation between the RPA variables and the single particle states. That is to say 

ocpj ocpj 

so that if one deals with collective electron motions, Eq. ^ can be approximated with very high degree of accuracy 
by the CP equations. 
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(56) 



This should be a reasonable approximation for many real molecular systems where one has delocalized valence electrons 
and relatively strong residual Coulomb interaction between them. We emphasize that our computational scheme is 
not restricted by the case of collective electronic excitation, i.e. by the approximation (Eq. ^). Taking into account 
the explicit expressionc3 for the functional derivatives of LOk in respect to the molecular orbitals 0i our dynamical 
equations can be straightforwardly used to do molecular dynamics on "non-collective" excited state surfaces. 

The dynamical equations for the RPA amplitudes are determined to be (assuming only a single RPA excitation in 
the system), 



{ph'\v\p'h)X^^,^,) 

(57) 



ph 



(5L^fc[0,X,y,R] 



p' h' 

+ ^kk^ph^ 



- E.uYX + E {{PP'\v\h'h)X%, - {ph'\v\p'h)Y^,,,) 



p' h' 



- Tfefel^t (58) 

or in matrix form 




A + r B \(\X) 

B A-rM|y) 



(59) 



One important point to consider is that the functional derivatives of Ehf with respect to the particle states (i.e. 
single particle states above the Fermi energy)vanish since the HF density matrix does not include these states. As we 
seen below, we need these states and their energies to compute the Coulomb matrix elements in the RPA dynamical 
equations. However, all is not lost since we have at hand the states below the Fermi energy and hence the ground 
state density. We can thus reconstruct and simply diagonalize Hi [p] to obtain the single particle excited states. 

Finally for the nuclear positions we have classical equations of motion on the excited state potential energy surface, 

5E,.,[c^,X,Y,Ii] 

5EHF[(t>,'R-n] 6uk[ct>,X,Y,Il^] 
(5R„ (5R„ 

The details of the derivation of each of the functional derivatives above are provided in the Appendix A. The derivation 
of the molecular dynamics equations of motion in the RPA compose the central result of this paper. 



IV. NUMERICAL AND ANALYTICAL EXAMPLES 



A. The dynamical optimization of the excited state on a two level model 



In order to get an indication of how our approach works, we apply our scheme to a model two-level system consisting 
of N distinguishable electrons labeled by p = 1,2, ...,7V. Each of these can occupy one of two orbitals, a lower and 
upper, having energies \£{R) and —^e{R) and distinguished by g = 1 and g = 2 respectively. The coupling constant 
V{R) does not depend on any quantum number. The gap between orbital e{K) and the strength of interaction V{R) 
depends on a classical variable R. The HF approximation is obtained by setting 1/ = 0, such that the lower level is 
fully occupied and the upper is empty, i.e. N ^ Q where i7 is the degeneracy of the levels. This model encapsulates 
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the salient-effects one expects to see in physically realistic systems. The model Hamiltonian possesses the SU(2) 
symmetrycJ and can be written in terms of quantum quasi-spin operators and the classical variable R: 



H = e{R)J^ - ^V{R) [j+J+ + J-J- 



(60) 



where the operators J+, J-, Jz are defined as follows 



J. 



p=i 



ip 



J_ = 



(61) 



Above a^pjCgp are particle creation and annihilation operators on the lower {g = 1) or the upper (g = 2) level. K is 
the spring constant which does not allow to classical system to move far from the equilibrium. 
Place the classical variable at the position 



R^Ro+q, 



(62) 



which is the sum of the equilibrium position Rq and the displacement q. Assuming the displacement from equilibrium 
is small, we can expand the Hamiltonian in q 



e{R)=eiR„)+^ 
m = ni?o)+^ 



q + Oiq^) 



Ro 



q + Oiq' 



(63) 
(64) 



Ro 



Truncating this at first order gives a linear interaction between classical and quantum systems. Neglecting the terms 
of 0{q^) we get the following quantum many-body Hamiltonian which depends parametrically on the displacement q 
of the classical variable: 



qj Jz 

q] (J+J+ + J-J- 



H(.)^f.(i.„)+^^^ 

Mq^ 1 2 
— — + -Kq^ 
2 2^ 



In the RPA, the operators j+, j_ behave like pure bosonic operators 



J_ J+ 

Jn' Vn 



= 1 



(65) 



(66) 



Using these bosons, we construct the RPA excitation creation operator (in the two-level model there is only one 
possible RPA excitation) 



(XJ+ -yj_) 



(67) 



The orthogonality of the RPA excited state imposes the following normalization condition on the RPA amplitudes: 
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(68) 



The RPA frequency is given by the matrix element 

uj[X,Y,q] = (RPA| [g, [Hiq),Q^] |RPA) 



e{Ro) 



de 
dR 



Ro 



q NXY 



(69) 



We can rewrite the u![X, Y, q] as a functional of only the Y amplitudes and classical q variable: 



.iY,q]=lsiRo)+^ 



g (i + 2y2) 



Ro 

dV 



N\/l + Y^Y 



(70) 



Ro 



- 2 \ViRo) 

The variational minimum of the RPA frequency can be determined analytically and is given as a the function of q: 

uj{q) 



\ 



Ro 



Ro 



The classical Lagrangian for the annealing of the electronic excited state is given by 

L[q, Y] = i/iy2 _ EHF[q] - c^[Y, q] + ^mq^ - ^Kq^ , 

where the EnFlq] is the HP ground state energy. 

The Lagrangian produces the equations of motion for the RPA amplitudes 



(71) 



(72) 



/^y = -4 eiRo) 



de 
dR 



q\Y 



2 V(«,.g 



q\N 



1 + 2Y^ 



Vl + Y^ 



mq 



de 
dR 



dV_ 
dR 



NVT+Y^Y - Kq 



(73) 



(74) 



The solution of these differential equations starting from an arbitrary value of Y{t = 0) and q{T — 0) is shown in 
Fig. 1 . The solutions oscillate about the minimum so that in order to get converged solution we introduced a velocity 
damping term to suppress oscillation around the equilibrium configuration. The final value of the amplitudes Y and 
X = vT+y^ are in excellent agreement with the analytical RPA results for a given equilibrium value of the classical 
displacement q^ ~ q{T oo): 



Y 



X = 



leiRo 


+ qo ) - 


^{qo) 


2uj{qa) 


le{Ro 


+ qo) + 


^{qo) 



2cj(qo) 



(75) 



Here the function e{Ro + qo) is given by Eq. |3| and uj{qo) by the Eq. 
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The simplified form of the usual CPMD equation of motion can be derived within our model if we set the residual 
interaction V{R) = and consider HF gap e(i?) as a fixed parameter. In this case only the classical degree of freedom 
evolve in time. One can see from the lower part of the Fig.l that there is a profound difference between excited and 
ground state dynamics for a classical system. This difference is due to the Hellmann-Feynman force from the RPA 
excited state (App. A.l). 



B. Excited state dynamics in conjugated polymer lattices 

In this example we take a more challenging case in which the ground state can be determined exactly, but the excited 
states cannot. The 7r-electrons in linear polyene systems, such as irans-poly-acetylene, can to a fair approximation 
be treated as a quasi-one dimensional electron gas within the tight-binding approximation. Accordingly, we start our 
consideration from the following model electronic Hamiltonian for the polymer chain 

H = Ho+: H,nt ■■ (76) 

Ho = + a("ri+l - Un))(cl+l,o-Cn.<T + c\ cjCn+l,a) 

+ fllK+i-"")' (77) 

11 

Hint ^\Y1 ™) X! cLcL^'Cmo-'Cna (78) 

nm aa' 

Here Cma creates an electrop of spin a on site m. The mean field Hamiltonian, Ho is based upon the so called 
Su-Schrieffer-Heeger (SSH)E3 model for the band-structure of trans-polyacetylene. The first term in the Ho gives the 
energy for an electron with spin a to hop between neighboring p-orbitals. The strength of this term is modulated 
by linear coupling to distortions in the polymer lattice away from evenly spaced lattice positions. Finally, the last 
terms in Ho gives the harmonic interactions between lattice-sites arising from the cr-bonds between neighboring lattice 
atoms. The : Hint '■ is the normal ordered, i.e. irreducible to one body operators, interaction between valence electrons 
from the p^-orbitals. 

The Hamiltonian Hg is quadratic form in fermion creation/annihilation operators. This quadratic form can be 
exactly diagonalized by the canonical transformation: 



^ k 



This canonical transformations is equivalent to a Hartree-Fock diagonalization of the electronic Hamiltonian and 
produces the uncorrelated mean field Hamiltonian: 

Ho = Y^ Ek{al„+aka+ - aL-Ofc,.-) + 2NKu' (79) 

where is the HF energy of a quasi-particle, aka±, with wavevector k in the first reduced Brillouin zone of the 
system, — 7r/2a < k < 7r/2a. We also define + to correspond to states above the Fermi-energy (particles) and — to 
those below (holes). The single-particle energy, Ek = y/f-l + A^, is a function of the unperturbed band energy, 
efc = 2tocos(fca), and the phonon- induced energy gap, Ak{u) = 2ausm{ka). A plot of the Hartree-Fock ground state 
energy surface in shown in Fig 2c. 

In this paper we are primarily interested in the details of computational algorithm and method. In order to make 
our consideration as model independent as possible we take the following general form of the residual interaction: 

H^nt = E E ^('?)«L±4'<T±afc-<z'T±afc'+9<T± , (80) 

kk' q aa' 

where the interaction V{q) is just a Fourier transformation of the physically meaningful interaction in coordinate 
representation (with U = 0.01 eV and Tq = a): 
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n,n+l 



1 + (a + 2M)/r, 



for even n 



n,n— 1 



1 + (a - 2u)/ro 
U for 71 even or odd 



for odd n, 



Vr. 



nn 



(81) 



Using the 2-body interaction of this form (with U = 0.01 eV and Tq = a) and the description of the ground state 
given above, we integrated-the RPA equations (Eq. ^ and Eq. ^) in /c-space over the first reduced Brillouin zone 
using the Verlet algorithmEJ for both the RPA amphtudes and classical coordinate. In all of our runs, we monitored 
the orthogonality of the amplitudes and re-orthogonalized whenever the orthogonality strayed from unity by 1 part 
in 10^. In the first example, we chose our initial excited state such that Xph corresponded to the first Cl(singles) 
triplet state and Yph = 0. In all of these calculations, the fictitious electronic mass was chosen to be /i=400 a.u. and 
the classical mass to be that of a -CH- monomer (m ~ 26000a.u.). 

In the first example, shown in Fig. 2a, we started from purely random distribution of initial amplitudes for X and 
Y again with u(0) = 0.1 A. As before, we dampened the RPA velocities to attempt to cool the system down to the 
excited state Born-Oppenheimer energy. Notice, that the system rapidly finds the lowest energy configuration of the 
excited state surface and begins to oscillate on the BO surface after only a few periods. This demonstrates the utility 
of our method in determining very rapidly minimal energy configurations on the excited state energy surface. 

In the next example, our initial excited state was taken to be an eigenstate of the A matrix (equivalent to the CI 
singles matrix) with orbital angular momentum 1 (triplet state). Starting from a lattice distortion of u=0.lA, we 
solved the full RPA equations of motion with the inclusion of a term to dampen the RPA amplitude velocities. This 
helped to keep the electronic velocities in sufficient check and prevented significant excursions from the excited state 
Born-Oppenheimer energy surface. The results of this calculation are shown in Fig. 2b. Here we have plotted the 
excitation energy, w as a function of u as the RPA amplitudes are quenched to their minimal values. Notice, however, 
that even though we kept the electronic velocities fairly "cool", there is some "fuzziness" to our energy surface. This 
is partially due to the fact that our initial state was not an exact solution to the RPA equations and to the fact that 
there is some dynamical coupling between the motions of the amplitudes and the classical degrees of freedom. Thus, 
the RPA amplitudes acquire a small amount of velocity as they are quenched to the final solutions. 

A final comment regarding the computational efficiency of our approach is in due order. First of all, in the 
calculations presented herein, we made use of an analytic representation of the ground state orbitals and energies. 
In general, one will need to determine these for each new molecular configuration either by using integrating the 
ground state CP equations or by converging the KS or HF equations. Assuming we were handed a set of single 
particle orbitals, we can compare the computational cost associated with each molecular dynamics step in the RPA 
with the cost required to setup and diagonalize the CI singles (CIS) matrix. In Fig. ^ we show the ratio of CPU 
times required to perform a CIS calculation versus the effort required to perform a similar sized single MD time-step 
using the RPA/MD method developed here as run on an SGI Origin2000 (on a single 195 MHz RIOOOO processor). In 
each case, the RPA/MD results were between 2 to 10 times as fast per time-step as diagonalizing the CI matrix with 
performance degrading with the size of the number of particle/hole states used. The CIS results reported are the CPU 
times required for the construction and diagonalization of the CIS matrix (converging only the lowest eigenstate) and 
does not include the additional effort need to compute analytical gradients which are needed in order to perform any 
sort of molecular dynamics on the excited state surface. Finally, we note that the each RPA/MD step requires the 
construct of both the A and B matrices whereas CIS requires only the A matrix. 

We also point out that the majority of the CPU effort per time-step went into the construction of the two-body 
interaction matrix elements. These matrix elements, along with the single particle orbitals, must be in hand at every 
time step and are a necessary component of any many-body treatment of the excited states. This effort scales as 
« 7V| due to the transformation between the quasi-particle basis and the lattice representation. For larger calculations 
involving the combination of many particle/hole state, this part of the calculation dominates the effort. Furthermore, 
in a full implementation of this scheme, one is required to diagonalize the HF or KS matrix to obtain the orbitals 
above the Fermi-energy. This too, adds to the net effort as it scales as « N^. In a general plane- wave implementation 
of this methods, we can utilize efficient and scalable multi-dimensional fast Fourier transformation algorithms and 
matrix-diagonalization methods to facilitate these two phases of the calculation. 



In this paper, we have presented a theoretical scheme for doing molecular dynamics simulations of electronically 
excited systems within the HF(KS)-|-random phase approximation. In the spirit of Car and Parrinello, we treat the 



V. DISCUSSION 
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RPA amplitudes, the single particle orbitals, and the nuclear coordinates as dynamical variables subject to a set of 
holonomic constraints which enforce the orthogonality relations of the RPA amplitudes and orbitals. 

As mentioned in the introduction, the two primary computational cost incurred in applying the CPMD scheme are 
enforcing the orthogonality constraints between the single particle states and the small time-step required to integrate 
the CP equations. As mentioned above, the constraint part of each iteration scales as ~ M x N'^. In the present 
scheme, we have and additional set of constraints for the orthogonality of the RPA states. Likewise, the evaluation 
of the RPA constraint should scale as ~ Mrpa x N'^pj^ where Mrpa is the number of particle-hole basis states used 
to construct the RPA excitation and Nrpa the number of RPA excitations in the system. Since we are primarily 
interested in systems with one RPA excitations and we can restrict the excitation to a limited number of particle- 
hole states, the added constraint does not pose a substantial increase in the computational effort. Furthermore, the 
orthogonality between the RPA state and the ground state is implicit in the construction of the RPA state and, hence, 
imposes no additional difhculties. Note, that this would not be the case had we chosen a perturbative construction 
of the excited state. 

A sample calculation on an albeit simple model system produces exact agreement with analytical results. While 
this model avoids a number of computational difficulties in using this approach, it does give an indication of how the 
equations of motion can be applied. The more realistic simulations of the excited state dynamics of one-dimension 
polymer chains demonstrate that our proposed approach is much more computationally advantageous since it con- 
siderably decreases the computational efforts over the direct solution of the RPA eigen-problem. We are currently 
developing a much more realistic and generalizable application of our excited states molepilar dynamics for simulating 
the photo-excitation dynamics in 7r-conjugated systems such as conducting polymers, Ea as well as the inclusion of 
transitions between RPA states, such as non-radiative triplet-singlet transitions. 
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APPENDIX A: EVALUATION OF FUNCTIONAL DERIVATIVES 

In this Appendix we compute the various functional derivatives of the RPA energy functional needed to solve the 
dynamical RPA equations presented above. 

1. Hellmann-Feynman forces for excited states 

We first invoke the Hellmann-Feynman theorerrSS to provide the forces for the evolution of the nuclei on the 
Born-Oppenheimer energy surface of the RPA excitations. 

The Hellmann-Feynman theorem is indeed valid for the RPA excited states, since 

<5(RPA|Q^,Qt jRPA) _ m^^\ [Q/.^,QIJ |RPA> _ 

SK - M - ° ^^^^ 

Thus, 

^i^..[0^X,r,R] ^ (RpAlQ^^^Qt^lRPA) (A3) 

To calculate this matrix element we first rewrite the gradient of H2 with respect to the nuclear coordinates in a second 
quantized form in the basis of electron creation and annihilation operators: 
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Here the matrix element 



(''m'^'^ = J '^^^^*WaR(^''"^'^^'^'"^ ^^^^ 

Using this, one can easily show 



(RPAIQ^,— Qt .|RPA) = (RPA| — |RPA) 



(RPA| 



|RPA) (A6) 



Neglecting the polarization fluctuation in HF vacuum, we arrive at the following expression for the classical forces: 

hh 

r) M 

php' 

phh' 

(A7) 

where the first term is simply the ground state Hellmann-Feynman force and the remaining terms are due to the 
electronic excitations. 

2. Functional derivatives for the amplitude motions. 

The functional derivatives for the evolution of the RPA state and HF vacuum are straightforward to compute and 
are presented below. 

a. The functional derivative of the HF energy with respect to the single particle/hole states. 

d2</.*(2)i/(21;p) 

y d40*(4))f](.^,|^|</',-) (A8) 



First, in vector form 



b. Functional derivatives for the RPA amplitudes 



^"^"^ A\Xk) + Bin) (A9) 



and 



6u>k 



S{Yk\ 

Written in terms of the components: 



B\Xk)+A\Yk). (AlO) 



15 



TTTfc; - / X^php'h'-^p'h' + ^php'h'Iph) (.Aiij 

°^ph p'h' 

'^^'^^t'yfcl ^ = X!("^P''P''»'-^p'''»' + Bphp'h'Xp^) (A12) 



APPENDIX B: DYNAMICAL CALCULATION OF THE LAGRANGE MULTIPLIERS FOR THE RPA 

EXCITED STATES 

The evolution of excited states is subject to holonomic constraints which enforce the orthogonahty of the RPA amph- 
tudes, ( Eq. ^9|). Because of this, we have constraint forces —TkkXpf^ and ^kkYph in the dynamical equations(Eq. |5^ 
and Eq. ^) which keep the RPA amplitudes on the constrained surface. However, the Lagrange multiplier F is 
unknown. 

Denote by 

ph 

Since ak is the integral of motion for the differential equations (Eq.^ and Eq. |5|) , the first and second time derivatives 
of ak should be zero. 

rj \ ^ vk^ yk xrk^\rk 
^ ^ph^ph ~ ^ph ^ph 

ph 

vk I x^fc* \^/c \/-k^\^k \^k^\^k \ {TiO\ 
^ph^ph ' ^ph^ph ~ -'ph -fp/i ~ -'ph -'p/i J \^^) 

We can then eliminate the "amplitudes accelerations" Y'^ ,X^, y^* ^X^* by the use of the RPA equation of motion. 
From these considerations, we obtain the expression for the Lagrange multiplier Tkk'- 

i fcfc — 




X 



o vk^ vk I \Ak^\rk 

^ l^ph ^ph^ph -Tp/i Jiph 

I^Ml^ph^ph '^ph'^ph) 



ph 



- Y — Y -4- V -4- V 



^Xl^ SX^h '"^SYX* 



(B3) 
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Figure Captions 



FIG. 1. The evolution of the RPA ampUtude Y (a.) and the displacement q (b.) in the classically-extended-SU(2) model. 
The dashed line is for simplified ground state CPMD; the result of our approach is shown by solid line. The following set of 
parameters is used in this calculation: fj, = M = l,e = 1, V = 0.06, TV = 10, {de/dr)Rg = 0.1, {dV/dR)Rg = 0.06, K = 10. 

FIG. 2. Energy surfaces for ground and excited states for linear polyene model, (a.) Excited state molecular 
dynamics starting from a random RPA vector converging to the excited state Born-Oppenheimer energy surface. The excited 

state energy surface is the sum of the excitation energy and the ground state surface plotted in c. Here, only the energy 
difference is plotted, (b.) Excited state molecular dynamics starting from the lowest triplet CI(S) eigenstate. As in a., only 
the energy difference between the triplet and ground state is shown, (c.) Ground state Born-Oppenheimer surface. Note the 
change in energy scale between the ground and excited state plots. See text for details. 

FIG. 3. Comparison of CPU effort for RPA/MD vs CI(S) per molecular dynamics time step £is a function of the number 
of particle/hole states. In both cases, the effort to construct the two body interactions is included; however, the CI(S) does 
not include the effort required to compute energy gradients. As the number of particle/hole states increases, this becomes the 

dominant part of the effort. 
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